!compute back projection imaging
module doa

   implicit none
   private

   public :: doaMUSIC
   public :: ptime
   public :: distance
   public :: spline1
!   public :: write_main

contains

!==================================================
! read in seismograms and parameters
!

subroutine doaMUSIC(ret)

use problem_class

type(problem_type),intent(inout) :: ret

integer :: ii,i,j,p,q,fli,fhi,lwork,INFO
double precision :: fc
complex(8), dimension(ret%n) :: a,ap,apUn
!complex(8), dimension(:),allocatable :: work
complex(8), dimension(10000) :: work1

complex*16, dimension(ret%n,ret%n) :: Rxx, unn, Un, Uv, Ae
!complex(8), dimension(ret%n,ret%ps*ret%qs) :: alib
real(8),dimension(ret%n,ret%n) :: Pm, Pm1
real(8), dimension(ret%n) :: W,ux,uy
real(8), dimension(10000) :: rwork

character(len=256) :: fileName

!allocate(a(ret%n))
!allocate(ap(ret%n))
!allocate(apUn(ret%n))
lwork=10000
!write(6,*)'work',work1(1:10)
!write(6,*)'rwork',rwork(1:10)

!allocate(rwork(lwork))
!allocate(work(lwork))
!do i=1,lwork
   !work(i)=(0,0)
   !rwork(i)=0d0
!enddo

!allocate(W(ret%n))
!allocate(Rxx(ret%n,ret%n))
!allocate(unn(ret%n,ret%n))
!allocate(Un(ret%n,ret%n))
!allocate(Uv(ret%n,ret%n))
!allocate(Ae(ret%n,ret%n))
!allocate(Pm(ret%ps,ret%qs))
!allocate(Pm1(ret%ps,ret%qs))
!allocate(Pw(ret%ps,ret%qs))
!allocate(ux(ret%ps))
do i=1,ret%ps
ux(i)=ret%uxRange(1)+ret%lat0+(ret%uxRange(2)-ret%uxRange(1))/(ret%ps-1)*(i-1)
enddo
!allocate(uy(ret%qs))
do i=1,ret%qs
uy(i)=ret%uyRange(1)+ret%lon0+(ret%uyRange(2)-ret%uyRange(1))/(ret%qs-1)*(i-1)
enddo

do i=1,ret%n
  do j=1,ret%n
       Pm(i,j)=0
  enddo
enddo
!fli=int(ret%fl/(ret%sr/ret%npts))+1;
!fhi=int(ret%fh/(ret%sr/ret%npts))+1;
fli=11
fhi=11
write(6,*)'fli=',fli
do ii=fli,fhi
  fc=(ii-1)*(real(ret%sr)/real(ret%npts))
  Rxx=ret%sxy(ii,:,:)
  ! write(6,*),'sxy',ret%sxy(:,1,2);

  Uv=Rxx
  
  write(6,*)'Rxx(2,:)',Uv(2,:)

   ! write(6,*)'W=',W
   write(6,*)'n=',ret%n
  do i=1,ret%n
    do j=1,i
      Uv(j,i)=conjg(Uv(i,j))
      if(i==j)then
      Uv(i,j)=real(Uv(i,j))
      endif
     enddo
  enddo
  
  !Uv(3,1)=3
  !Uv(1,3)=3
  !write(6,*)'ori',Uv
  !do i=1,lwork
  ! write(6,*)'i',i
  ! work1(i)=(0.0,0.0)
  ! rwork(i)=0d0
 ! enddo

  call zheev('V','U',ret%n,Uv,ret%n,W,work1,lwork,rwork,INFO)
  write(6,*)'eigs=',W(ret%n-5:ret%n)
  write(6,*)'Uv',Uv(:,1)
  do i=1,ret%n
    do j=1,ret%n
       if(i==j)then
       Ae(i,j)=W(i)
       else
       Ae(i,j)=(0.0,0.0)
       endif
       unn(i,j)=(0.0,0.0)
       Pm1(i,j)=0.0
     enddo
   enddo
   do i=1,ret%n
      do j=1,ret%n-ret%Mm
         unn(i,j)=Uv(i,j)
      enddo
   enddo
   write(6,*)'n-m',ret%n-ret%Mm
     write(6,*)'n',ret%n
   write(6,*)'Mm',ret%Mm

   write(6,*)'Uv',Uv(1:2,1:2)
   write(6,*)'unn',unn(1,1:10)
   Un=matmul(conjg(unn),transpose(unn))
   write(6,*)'Un',Un(1:2,1:2)
   write(6,*)'info',INFO
   do p=1,ret%ps
     do q=1,ret%qs
       call ptime(ret%n,a,fc,ret%lat0,ret%lon0,ux(p),uy(q),ret%lat,ret%lon,ret%rtn,ret%rr,ret%tt)
      ret%alib(:,p,q)=a
      ! write(6,*)'a=',a(1:10)
      ! write(6,*)'lat0,lo0,uxp,uyq,uxrange,uyrange',ret%lat0,ret%lon0,ux(p),uy(q),ret%uxRange,ret%uyRange
     enddo
   enddo

   do p=1,ret%ps
     do q=1,ret%qs
      a=ret%alib(:,p,q)
       ap=conjg(a)
       !do j=1,ret%n
       !apUn(j)=dot_product(ap,Un(:,j))
       !enddo
       apUn=matmul(a,Un)
      !apUn=matmul(a,Rxx)
      !write(6,*)'Un',Un(1,1:10)
       Pm1(p,q)=abs(dot_product(a,a)/dot_product(apUn,a))
      ! Pm1(p,q)=abs(dot_product(apUn,a)/(dot_product(a,a)))
       ! Pm1(p,q)=abs(dot_product(a,a))
     enddo
   enddo
   write(6,*)'Pm',Pm1(1,1)
   Pm=Pm+Pm1
enddo
write(6,*)'ok before write file'
write(6,*)'Pm',Pm(1,1)
write(fileName,'(f0.1)')ret%t1
fileName='Pm'//fileName
!write(*,*)'Pm',Trim(fileName)
!fileName='Pm':q

write(6,*)fileName
open(unit=121, FILE=fileName)

do p=1,ret%ps
 do q=1,ret%qs
!  write(6,*)Pm(p,q)
  write(121,'(f0.5a)',advance='no')Pm(p,q),' '
!  wirte(17,*)' '
 enddo
  write(121,"(a)")' '
enddo



write(6,*)'ok after write file'

end subroutine doaMUSIC



!============================= ptime compute the travel time difference 
subroutine ptime(n,a,fc,lat0,lon0,uxp,uyq,lat,lon,rtn,rr,tt)

integer, intent(in) ::n,rtn
double precision, intent(in) ::fc,lat0,lon0,uxp,uyq
double precision, dimension(n), intent(in) :: lat,lon
double precision, dimension(rtn), intent(in) :: rr,tt
complex(8), dimension(n), intent(out) :: a
integer :: i
double precision, dimension(n) :: rtmp,r1, r2, t1, t2, t
do i=1,n
call distance(r1(i),lat0,lon0,lat(i),lon(i))
call distance(r2(i),uxp,uyq,lat(i),lon(i))
call spline1(rtn,rr,tt,r1(i),t1(i))
call spline1(rtn,rr,tt,r2(i),t2(i))
enddo

t=t2-t1;
t=t-sum(t)/real(n);
!write(6,*)'t=',t(1:10)
!write(6,*)'fc',fc
do i=1,n
a(i)=exp(-(0.0,1.0)*2*3.1415926*fc*t(i));
enddo

end subroutine ptime
!=========================spline1 does 1d spline interploation
subroutine spline1(n,t,y,x,s)
      integer, intent(in) :: n
      double precision, dimension (n) :: t, y
      double precision, intent(in):: x
      double precision, intent(out):: s     
      integer :: i
      double precision :: diff
      do i = n-1,1,-1     
       diff = x - t(i)     
        if (diff >= 0.0) exit  
      end do    
      s = y(i) + (x - t(i))*(y(i+1) - y(i))/(t(i+1) - t(i))
end subroutine spline1
!========================distance computes the great circle distance in degrees
subroutine distance(d,lat1,lon1,lat2,lon2)

      double precision, intent(in):: lat1,lon1,lat2,lon2
      double precision :: a
      double precision, intent(out) :: d
      a=dsind((lat2-lat1)/2)**2+dcosd(lat1)*dcosd(lat2)*dsind((lon2-lon1)/2)**2
      d=180/3.1415926*2*atan2(sqrt(a),sqrt(1-a))
end subroutine

end module doa
